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PREFACE 


This  report  describes  the  development  and  application  of  a  numerical 
transport  model  used  as  a  basis  for  2-D  vertically  averaged  estuarine  water 
quality  models.  The  preparation  of  this  report  was  sponsored  by  the  Office, 
Chief  of  Engineers  (OCE),  under  the  Environmental  Impact  Research  Program 
(EIRP).  The  Mobile  District,  CE,  sponsored  the  Mississippi  Sound  Study,  which, 
is  presented  as  a  test  application.  Technical  Monitors  for  EIRP  were  Dr.  John 
Bushman  and  Mr.  Earl  Eiker  of  OCE  and  Mr.  David  B.  Mathis,  U.  S.  Army  Water 
Resources  Support  Center. 

The  work  presented  in  the  report  was  conducted  from  July  1979  through 
June  1983  in  the  Wave  Dynamics  Division  (WDD)  of  the  Hydraulics  Laboratory  (HL) 
of  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES)  under  the  general 
supervision  of  Messrs.  H.  B.  Simmons,  Chief,  HL,  and  Claude  E.  Chatham,  Jr., 
Acting  Chief,  WDD.  The  WDD  was  transferred  to  the  Coastal  Engineering  Research 
Center  (CERC)  of  WES  on  1  July  1983.  From  July  through  September  1983,  work 
was  performed  in  the  WDD  of  CERC  under  the  general  supervision  of  Dr.  R.  W. 
Whalin,  Chief,  CERC,  and  Mr.  Chatham,  Chief,  WDD.  Dr.  R.  A.  Schmalz,  Jr.,  WDD, 
conducted  the  Mississippi  Sound  Study  and  prepared  this  report. 

The  preparation  of  this  report  was  monitored  by  Mr.  Ross  W.  Hall,  Ecosys¬ 
tem  Research  and  Simulation  Division  (ERSD) ,  Environmental  Laboratory  (EL), 
under  the  general  supervision  of  Mr.  Don  L.  Robey,  Chief,  ERSD,  and  Dr.  John 
Harrison,  Chief,  EL.  Program  Manager  at  WES  for  EIRP  was  Dr.  Roger  T.  Saucier. 

Commanders  and  Directors  of  WES  during  the  conduct  of  this  study  and  the 
preparation  and  publication  of  this  report  were  COL  Nelson  P.  Conover,  CE,  and 
COL  Tilford  C.  Creel,  CE.  Technical  Director  was  Mr.  F.  R.  Brown. 
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Two-Dimensional  Vertically  Integrated,  Time-Varying  Estuarine 

Transport  Model,"  Instruction  Report  EL-85-1,  US  Army 

Engineer  Waterways  Experiment  Station,  Vicksburg,  Miss. 
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CONVERSION  FACTORS,  U.  S.  CUSTOMARY  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


t,  S 

met  r 

t  eet 
tfft 
m .  I  f 


.  customary  units  of  measurement  used  in  this  report  can  be  converted  to 
ir  (SJ)  units  as  follows: 


Multiply 


per  second 

s  per  hour  (U.  S.  statute) 


By 

0.3048 

0.3048 

1.609347 


To  Obtain 
metres 

metres  per  second 
kilometres  per  hour 


USER  GUIDE  FOR  WIFM-SAL:  A  TWO-DIMENSIONAL  VERTICALLY 
INTEGRATED,  TIME-VARYING  ESTUARINE  TRANSPORT  MODEL 

PART  I:  CAPABILITIES  AND  LIMITATIONS 

1.  The  transport  model  WIFM-SAL  was  developed  as  a  prerequisite  require¬ 
ment  for  water  quality  models  to  be  used  in  the  analysis  of  water  quality  prob¬ 
lems  in  shallow  estuaries  and  embayments  which  may  be  considered  vertically 
well  mixed  thereby  justifying  a  vertically  integrated  approach.  The  model  is 
two-dimensional  in  the  horizontal  and  generates  time-varying  water  surface 
elevations,  velocities,  and  constituent  fields  over  a  space  staggered  grid. 
Units  of  measure  are  expressed  in  the  English  system  (slug-foot-second). 

2.  Two  constituent  transport  schemes  have  been  incorporated  in  the 

U.  S.  Army  Engineer  Waterways  Experiment  Station  (WES)  Implicit  Flooding  Model 
(WIFM)  developed  by  Butler  (1980).  Constituent  computations  are  performed  at 
the  same  time  step  interval  as  employed  in  the  hydrodynamic  computations. 
Therefore,  if  desired,  the  user  may  develop  the  coding  necessary  to  density 
couple  the  hydrodynamics  if  this  is  important  for  the  problem  of  concern. 
Density  coupling  is  not  implemented  in  the  model  at  this  time. 

3.  An  exponentially  stretched  grid  system  is  used  in  WIFM-SAL  allowing 
the  user  to  increase  resolution  in  specific  areas  where  more  computational  de¬ 
tail  is  desired.  This  feature  is  particularly  useful  in  modeling  inlets  and 
barrier  island  systems. 

A.  Since  the  constituent  transport  schemes  are  directly  encoded  within 
WIFM,  this  model  must  be  used  to  provide  the  hydrodynamic  description.  Future 
work  will  be  conducted  to  develop  a  separate  transport-dispersion  model,  allow¬ 
ing  for  user  selectable  hydrodynamic  input  and  transport  scheme  selection. 

5.  Although  WIFM  has  been  used  extensively  in  moving  boundary  applica¬ 
tions,  the  transport  schemes  assume  a  fixed  land/sea  boundary.  Future  work  is 
needed  to  remove  this  restriction. 

6.  The  constituent  transport  equation  considered  is  for  a  passive 
scalar  without  source/sink  terms.  The  extension  to  multiple  (reacting)  con¬ 
stituent  systems  remains  to  be  developed. 

7.  WIFM-SAL  allows  the  model  user  to  employ  results  computed  on  a  global 
grid  as  boundary  conditions  on  a  more  spatially  limited,  refined  grid 


concentrated  around  the  area  of  interest.  In  addition,  the  user  may  select 
either  of  two  distinct  transport  schemes.  Scheme  1  is  a  flux-corrected  trans¬ 
port  scheme  capable  of  resolving  sharp  fronts  without  oscillation.  Scheme  2 
is  a  full,  three  time  level  scheme  directly  compatible  with  the  three  time 
level  hydrodynamics.  Scheme  1  requires  approximately  three  times  more  com¬ 
puter  time  than  scheme  2  but  is  more  accurate  than  scheme  2  for  sharp  front 
problems. 

8.  However,  on  a  coarse  spatial  resolution  global  grid  covering  a  large 
area,  scheme  2  results  may  be  used  in  areas  away  from  sharp  fronts  to  provide 
boundary  conditions  for  a  more  refined  grid  system  encompassing  the  sharp 
front  region  of  propagation.  Scheme  1  may  then  be  selected  to  resolve  the 
sharp  front  over  this  refined  grid.  Thus,  the  telescoping  grid  capability  in 
conjunction  with  the  user  selectable  constituent  transport  scheme  is  an  ex¬ 
tremely  powerful  concept  in  real  world  transport  problem  solving. 


PART  IT:  THEORETICAL  DEVELOPMENTS 


9.  Consider  the  instantaneous  three-dimensional  constituent  transport 
equation.  The  time  scales  for  which  this  equation  applies  are  of  much  shorter 
duration  than  can  be  modeled.  Therefore,  the  instantaneous  equation  is  tem¬ 
porally  averaged.  Under  the  vertically  integrated  approach,  the  resulting 
equation  is  then  depth  averaged.  The  transport  equation  obtained  is  then 
transformed  using  an  exponential  stretch.  Numerical  approximations  to  the 
transformed  equation  are  formulated  followed  by  the  development  of  relations 
for  the  effective  dispersion  coefficients. 

Constituent  Transport  Equation  in  Cartesian  Coordinates 

10.  The  instantaneous  constituent  transport  equation  is 

3s  3s  .  3s  3s  3  /  3s  \ 

3t  U  3x  V  3y  W  3z  3x  \  x  3x  ) 

♦fe  MW;  k  I) 


x,y,z  =  Cartesian  coordinates 

u,v,w  a  velocity  components  in  the  x-  ,  y-  ,  and  z-directions , 
respectively 

t  =  time 

s  =  concentration  of  the  material  of  concern 

D  s  molecular  diffusion  coefficient  in  the  x-direction 
x 

D  =  molecular  diffusion  coefficient  in  the  y-direction 

y 

D  =  molecular  diffusion  coefficient  in  the  z-direction 
z 

For  a  turbulent  flow,  the  turbulent  diffusion  is  much  greater  than  the  molecu¬ 
lar  diffusion.  The  following  analogous  formula  holds  where  time  averaging 
over  the  time  scale  of  the  turbulence  has  been  performed. 


3s  3s  ,  3s  3s 

-  +  11  —  +  v  —  +  w  • — 

3t  3x  3y  3z 


(K*l) 


9  1 

rK  9s\ 

+  ( 

3y  ' 

vKy  9y/ 
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\z  3z/ 

where  K  ,  K  ,  and  K  are  turbulent  diffusion  coefficients.  Equation  2 
x  y  z  ^ 

may  be  written  in  conservation  form  by  adding  s  times  the  continuity  equa¬ 
tion  (namely,  zero)  to  the  left-hand  side  to  obtain 


9s  +  3(us)  +  9 (vs)  9(ws) 

9t  9x  9y  3z 


9_ 

3x 


(3) 


This  form  of  the  equation  is  then  depth  integrated  as  described  in  Schmalz 
(1981a)  to  obtain: 


9_ 

9t 

where  h  is 
coefficients 


(hs)  +  |j  (hus)  *  ly  (hvs)  =  §j  (hK*  |l)  ♦  |y  (hK*  I) 
the  water  depth  and  K*  and  K*  are  effective  dispersion 


(A) 


Constituent  Transport  Equation  in  Transformed  Coordinates 


11.  The  transport  equation  is  transformed  from  x  -  y  space  to 
-  (*2  space  by  means  of  the  following  coordinate  transformation  as  con¬ 
sidered  by  Butler  (1980). 


x 


a!  +  Vi 


(5) 


y  =  a2  +  b2a2 


(6) 


The  terms  a|  *  b  ,  c^  ,  a2  ’  ^2  ’  an<*  C2  are  constants  valid  for  dif¬ 
ferent  regions  in  the  grid.  Then  for  an  arbitrary  hydrodynamic  variable 
P(x,y,t) 


9p  _  9p  dai 
9x  9a  dx 


9p  _  9p  da2 
3y  9^2  dy 
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If  we  introduce  px  =  dx/dc^  and  p2  =  dy/dc<2  then 

3£  _  1_  3g_  §£  _  1_  3£_ 
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Considering  Equation  8a  in  an  alternate  manner 
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Employing  previous  notation,  Equation  13  is  rewritten  as  follows 
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we  obtain 


Note,  however,  from  the  relation  between  3/3x  and  d/da^ 

*  S^e/L.)2  ,  ?e_  _d_ /j_\  i_ 
a*2  a^VV  3“i  d“iVVtJi 

This  relation  is  equivalent  to  Equation  S. 

12.  If  we  consider  a  hydrodynamic  variable  p(aj,a2,t)  and  let 
j*  ,  n  be  defined  such  that 


Pi*,j*  =  P(i'v^a2,j*A0'l’nAt) 

Then  let  i  ,  j  ,  n  be  such  that 

Pi,j  =  P[a2  +  b2(i*Aa2)  2  *  al  +  b!  C  j*Aa1 )  1  ,  nAtj 

We  employ  uniform  spacing  in  a^  -  a space  and  irregular  .spacing  in 
space.  We  may  evaluate  the  derivatives  with  respect  to  x  and  y  as 
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where 
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For  the  second  derivative  term  we  obtain 
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where 
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=  h(j*Aa1 ) 


Similarly,  for  — ^ 

dy 


The  underlined  terms  in  Equations  10  and  11, 


i»J 


although  they  may  be  computed  exactly,  are  approximated  using  finite  differ¬ 
encing  on  and  P2  . 

13.  Transforming  Equation  4  in  x  -  y  space  to  Oj  -  a 2  space  we  ob¬ 
tain  the  following  result: 
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where  d  is  introduced  as  the  depth  in  place  of  h 

(  )t  =  a/at 
(  )„  =  a/aa. 


a. 


1 


(  )fl2  =  3/3«2 

Equation  21  is  the  relation  that  is  the  subject  of  numerical  approximation. 


Numerical  Approximations 


14.  Schmalz  (1983a,  1983b,  1983c)  considered  several  alternate  tech¬ 
niques  for  approximating  Equation  21.  The  Flux  Corrected  Transport  Scheme 
(FCT)  was  selected  as  the  most  accurate  scheme  and  has  been  incorporated  in 
the  Waterways  Experiment  Station  Implicit  Flooding  Model  (WIFM) .  In  addition 
a  three  time  level  explicit  transport  scheme  was  also  incorporated  in  the 
model.  A  space  staggered  grid  as  shown  in  Figure  1  was  employed  in  all  of  the 
formulations.  The  datum  convention  is  presented  in  Figure  2. 

15.  Let  us  introduce  the  following  notation  as  a  prelude  to  the  approxi- 

mations.  Define  for  an  arbitrary  variable  F  ,  where  t  =  kAt  ,  y  =  nAy 

n  ,m 

x  =  mAx  : 


ffk  )  = 

:  pk+1/2  _  fk 

(22a) 

\  n,m/ 

n,m  n,m 

) 

=  Fk+1  -  Fk 

(22b) 

t\  n,m/ 

n,m  n  ,m 

(pk  ) 
j \  n,m/ 

=  Fk  -  Fk 

n,m+l/2  n,m-l/2 

(22c) 

5  (Fk  )  = 

a  \  n,m/ 


F^  -  Fk 

n+l/2,m  n-l/2,m 


a  /pk  Fk  \ 

~  .  Vn.m+1/2  n,m-l/2/ 


n,m 


h  (f“  .  +  F 


rk 

n+l/2,m 


K 

n-l/2,m/ 


(22d) 


(22e) 


(nor  \ 


- * - 

16.  Two  schemes  are  used  in  implementing  this  approach:  a  lower  order 
in  space  nonoscillatory  scheme  and  a  higher  order  in  space  scheme  subject  to 
oscillation.  In  the  method  implemented,  two  time  level  implicit  multiopera- 
tional  ADI  schemes  were  employed.  The  forward  time  upwind  space  (FTUS)  and 
forward  time  centered  space  (FTCS)  schemes  were  used  as  the  lower  and  higher 
order  in  space  schemes,  respectively,  and  are  discussed  in  turn  below.  Fi¬ 
nally,  the  necessary  flux  correction  procedures  are  developed. 

17.  Leendertse  FTCS  multioperational  scheme.  The  following  finite  dif¬ 
ference  equation  is  considered  as  an  approximation  to  the  nonlinear  transport 
equation  (21): 


6 '  k(ds)  +  — — 7—  6 

t  2Ao1 (p. )  a. 

1  i  m  I 


(a  a  a  a  \ 

V— jk+1— k+1  k+1  ^  — jk— k  kJ 
\d  s  u  +  dsu/ 


.  At  * 
2Aa2(p2)n  % 


(a_  a„  a _  a9  \ 

-k+1— k+1  k+1  ^  -jk— k  k  ) 
d  s  v  +  d  s  v  / 


a.  6  (sk+1)  a  6  (sk) 

At  .  -k+1  .k+1  “l  ~~ k  k  “l 

2(Ao  )2(y  )  “l  “l  “l  (yl}m 

JL  1  m  L  J 

a.  6  (sk+1)  a?  6  (sk) 

At  .  -k+i  „k+i  a2  -k  „k  a2 

- -  6  a  K.  7  r  +  d  K  \ 

2(Aa„)2(y9)  a2  “2  (u2)n  a2  (u2)n 

l  L  n  *- 


=  0  at  (n,m) 


The  solution  of  the  above  semi-implicit,  difference  scheme  requires  the  inver¬ 
sion  of  a  large  unbanded  matrix.  In  order  to  reduce  computational  effort,  the 
following  ADI  multioperational  difference  equations  are  used. 

18.  The  approximations  for  the  X-Sweep  may  now  be  written  as  follows: 


At  5  /  “l  ttl  \ 

al  \  .k+1/2*  k+1/2*  k+1/2*/ 

+  o«."  77. . r  'd  s  u  / 


^<“S)  '  2tol(„,). 
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x  *„  Csk+1/2*) 
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a  k. 


a 


(y  i ) 
l  m 


(2  i 


a2  a2 
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2(p2)nAa2  a2 


,k  k  k 
yd  s  v 


At  6 


ai 


2Aa2(n2)n 


“i 

-jk  k  a2 

a  K 


6  (sk) 


a2  (Vn 


=  0  at  (n,m) 


If  we  place  all  terms  at  time  level  k+1/2*  on  the  left-hand  side  of  the 


equation  and  expand  =  Kq 


, ...k+i/2*  at 

n.m  Zaa^Uj^ 


V  k+1/2*  _ 

\nntm»l  n.m+1 


'i 

k+1/2*  \  /  k+1/2*  k+l/2*\ 
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24oi(ui)« 
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■  <n~l _ n,ni-l  nn.m _ n.tn/  V  sn,m _ 8n,a-l  /  ^k+1/2* 

2 -  (ul,a-l/2  *n,a-l/2. 


Collecting  all  terms  in  Equation  23  at  time  level  k  denoting  the  result  as 

B  ,  we  obtain  K  =  K 
m  y  a2 


”  «  *  ,•*•»>  ,  »  »  «  a  .  «  .  •  L  *  W  •  ^  *  *  *  .  •  «  *  .^  I 

■■*--  *--*-*  At a  ■»  «.  *-  «-  ■ 
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B  -  (ds) 


/  k  k  \ 

/  k  k  \ 

At 

\nn+l,m  n+l.m  nn,m  n.m/  k 
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\  n+l,m  n.m/ 
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\  n-l.m  n-l»m  n.m  n.m/  k  \  g-l,m  n,m/ 

2.  n-l/2,m  2. 
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(26) 


B 


In  Equation  25  we  define  -a  ,  ,  a  , ,  ,  and  a  as  follows 

n,m-l  n,m+l  n,m 


/a. 


-a 


n,m-l 


AtH^* 

_ n,m-l/. 

2Aa. (u  ) 
i  i  m 
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(28) 


At 


a  =  dk+l/2*  + 

n,m  n,ra  2Aa, (p, ) 

l  i  m 


[^k+1iu2 

n , m+1 / 2 
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At 
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(29) 


19.  Collecting  all  results  we  obtain  the  following  interior  equa¬ 
tion  for  the  X-Sweep 


k+1/2*  k+1/2*  A 

a  ,  s  '  +  a  s  + 

n,m- l  n,m- I  n,m  n,m 


gk+l/2*  =  B 
n,m+l  n,m+l  m 


(30) 
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20.  The  approximations  for  the  Y-Swcep  may  now  be  written  as  follows: 


At6  (  2  2  \ 

6k+1/2*(ds)  ♦  2  \Hk+1  «k+1  k+1) 
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At6 


1 
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1  m 


(31 


=  0  at  (n,ra) 


Expanding  Equation  31  by  employing  Equation  22  and  collecting  terms  at  time 
level  k+1  on  the  left-hand  side  and  leaving  terms  at  time  level  k+1/2*  on 
the  right-hand  side,  the  following  interior  equation  for  the  Y-Sweep  is 
obtained: 


a  sk+1  +  a  sk+1  +  k+1 

n-l,m  n-l,m  n,m  Sn,m  an+l,m  Sn+l,m  ~  n 
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21.  Leendertse  FTUS  multioperational  scheme.  The  following  finite  dif 
ference  equation  is  considered  as  an  approximation  to  the  nonlinear  transport 
equation  (21): 


/  a1  k+1  a1  k 
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22.  The  following  upwind  difference  operators  are  used  in  the 
above  equation  and  are  defined  at  (n,m)  as  follows: 
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*n ,m-l/2 
k 

*n ,m+l /2 
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23.  To  effect  the  solution  of  this  scheme,  the  inversion  of  an  unbanded 
matrix  is  again  required.  Thus,  an  ADI  scheme  similar  to  the  previous  tech¬ 
nique  (upwind  differencing  is  employed  for  the  advective  terms)  is  used.  The 
necessary  modifications  for  the  X-Sweep  are  shown  in  Table  1  while  those  em¬ 
ployed  for  the  Y-Sweep  are  given  in  Table  2. 

24.  Flux  correction  procedures.  If  the  factorization  terms  are  ignored, 
the  schemes  above  may  be  written  in  the  following  flux  format: 


dk+1sI  =  dk  sk  -  ["aci  (p  )  Aot.(y-)  ~|  ^(f1 
n,m  n,m  n,m  n,m  L  1  1  m  2  2  nj  \n 


n+1/2  ^jn 


-  F1  +  F1  -  F1 

n-l/2,m  n,rrrt-l/2  n,m-l/2> 


where  t  =  kAt  ,  x  =  £  (lO.Act.  ,  y  =£  (y2).Aa„ 

i  x  1  i 

S  =  concentration  at  location  (n,m)  at  time  level  k 
n,m 

Aa^u,)  =  x  space  step  at  m 
11m 

^°2^2^n  E  y  space  step  at  n 

I  =  general  index  at  time  level  k+1  ,  which  we  set  to  H  or 
L  for  the  higher  or  lower  scheme,  respectively 

j 

F  +  ._  +i/2  “  fluxes  through  the  appropriate  cell  faces  of  cell  (n,m). 

'  ’  ~  '  Form  is  dependent  upon  the  finite  difference  formulation 

We  observe  from  Equation  39  that  the  difference  between  the  higher  and  lower 
order  scheme  at  (n,m)  may  be  written  as  follows: 

(Sn,m  "  Sn,m)“  [(Fn+l/2,m  '  Fn+l/2,m) 

-  (Fn-l/2,ra  -  Fn-l/2,m)  +  (Fn, n+1/2  "  Fn,n+l/2)  (40) 


Fn, m-l/2)J 


Note  this  difference  may  be  expressed  as  an  array  of  fluxes  between  adjacent 
grid  points  and  is  the  condition  required  to  effect  the  flux  correction  pro¬ 
cedures  as  given  by  Zalesak  (1979).  We  next  develop  the  flux  expressions  for 
the  higher  (F  )  and  lower  (F  )  order  schemes.  In  order  to  aid  in  notation, 
we  make  the  following  definition  for  an  arbitrary  variable,  F  : 


Table  1 

X-Sweep  Modifications  FTUS 


Equation  FTCS  FTUS 


25.  For  the  higher  order  scheme  we  employ  the  FTCS  scheme  written  in 
Equation  23  in  which  the  factorization  terms  developed  in  the  multioperational 
method  are  not  shown.  Equation  23  may  be  written  in  the  form  of  Equation  39, 
where  the  total  fluxes  are  presented  as  the  sum  of  advective  and  diffusive 
fluxes . 

26.  From  Equation  23  one  then  obtains  for  the  advective  fluxes: 
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(42) 


(43) 


The  diffusive  fluxes  are  then  given  by  the  following  relations 

(K  =  K  ,  K  h  K  ): 
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27 .  For  the  lower  order  scheme,  the  FTUS  scheme  written  in  Equation  37 
is  employed.  Factorization  terms  generated  by  the  multioperational  method  are 
not  considered.  Equation  37  is  written  in  the  form  of  Equation  39.  The  total 
fluxes  are  presented  as  the  sum  of  advective  and  diffusive  fluxes. 

28.  Fronr.  Equation  37  one  obtains  the  following  set  of  advective  fluxes: 
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The  diffusive  fluxes  are  obtained  from  Equations  44  and  45  with  H  replaced 
by  L  . 

29.  The  antidif fusive  fluxes  are  then  computed  as  follows: 
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In  computing  the  difference  between  the  diffusive  fluxes  (third  and  fourth 

terms  in  the  above  expressions),  note  that  the  terms  with  Sk  ma^  COm 

n,m 

pletely  eliminated. 

30.  Next  the  maximum  and  minimum  cell  values  are  determined: 
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31.  Next  the  sum  of  all  antidif fusive  fluxes  into  cell  (n,m),  P  , 

n.m 


is  determined: 
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The  maximum  allowable  mass  into  the  cell,  Q  ,  is  then  computed  as  follows: 
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(56 


32.  Similarly,  the  sum  of  all  antidif fusive  fluxes  out  of  cell  (n,m), 
is  determined: 


P 

n,m 


P  =  max(0,A  , .  )  -  min(0,A  n  ) 

n,m  *  n+l/2,m  n-l/2,m 


+  max(0‘An,m+l/2)  "  min(°’An,m-l/2)  (57 


The  maximum  allowable  mass  to  leave  the  cell,  Q  ,  is  then  computed: 

n,m 


Q"  =  (sL  -  Smin){(M  )  Aa  (|J  )  Aa  dk+1l 
n,m  \  n,m  n,m/  [  1  m  1  M2  n  2  n,mj 


33.  The  following  ratios  are  next  computed  for  use  in  determining  the 
limiting  coefficients: 


(l,Q+  /  P+  )  P+  >  0 

\  n,m  /  n,m  /  n,m 


P+  =  0 
n,m 


[min  (l,Q~  /  P‘  J  P”  >  0 

I  \  n,m  /  n,m/  n,m 


P  =  0 
n,m 


The  limiting  coefficients  are  then  given  by 


'n+1/2 ,m 


(min  (r+..  ,R  )  A  >  O] 

J  \  n+l,m  n,m/  n+l/2,m  -  I 

/  min  (r+  ,R  ..  )  A  m  <  0i 

\  \  n,m’  n+l,m/  n+l/2,m  ) 


"n , m+ 1 / 2 


lin  (r+  ..,R"  \  A  >  o) 

\  n,m+l  n,m/  n,m+l/2  -  { 

lin  (r+  ,R  ..)  A  <  Oi 

\  n,m  n,m+l/  n,m+l/2  / 


34.  The  antidif fusive  fluxes  in  Equations  50  and  51  are  limited  by 
multiplying  by  the  limiting  coefficients  and  the  solution  is  advanced  to  the 
next  time  level: 


Sk+1  =  SL  -  [aot  (p  )  Aa  (|J  )  dk+1l  C 
n,m  n,m  |_  11m  2  2  n  n,nu 


n,m  n,m  [_  11m  2  2  n  n,nu  n+l/2,m  n+l/2,m 

^n-l/2,m^n-l/2,m  +  <"n,m+l/2^n,m+l/2  ^n,m-l/2^n,m-l/2  (62) 

k+1  L 

We  observe  that  for  C  . ,  =  C  , ,  =  0  ,  S  =  S  and  for 

n+l/2,m  n,m±l/2  n,m  n,m 

C  =  C  =10  Sk+^  = 

n±l/2,m  n,m±l/2  •  n,m  n,m  ' 

35.  The  coding  of  the  flux  corrected  transport  procedures  is  presented 
in  Subroutine  CONC  in  Appendix  A. 

Three  time  level  explicit  scheme 

36.  In  order  to  avoid  the  averaging  of  hydrodynamic  quantities,  which 
is  performed  when  employing  a  two  time  level  transport  scheme  with  a  three 
time  level  velocity  scheme,  a  three  time  level  explicit  scheme  is  considered. 

37.  It  is  instructive  to  observe  the  form  of  the  continuity  equation 
employed  in  the  multioperational  hydrodynamic  scheme. 

X-Sweep : 
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(64  j 


at  (n,m) 


where 


At  =  time  step  length 


H*  =  water  surface  elevation  at  intermediate  time  level  * 
k±l 

n  =  water  surface  elevation  at  time  level  k±l  at  cell  (n,m) 
n,m  ’ 


AOj  =  space  increment 

Aoi^  =  space  increment 


k+1 


un  m+1/2  =  x  “  velocity  component  at  time  level  k+1  at  cell  (n,m) 


k-1 


un  m+1/2  =  x  "  velocity  component  at  time  level  k-1  at  cell  (n,m) 
k+1 


vn+l/2  m  =  y  ’  a2  component  at  time  level  k+1  at  cell  (n,m) 

k-1 


/n+l/2  m  5  y  ’  a2  vel°city  component  at  time  level  k-1  at  cell  (n,m) 


d  =  water  depth  at  time  level  k  at  cell  (n,m) 
n  jin 


If  we  eliminate  the  intermediate  level  r|*  ;  e.g.,  solve  for  r|*  in  Equa¬ 
tion  63  and  substitute  in  Equation  64,  we  obtain: 
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at  (n,n) 


k  k 

Since  d  =  n  -  h  ,  Equation  65  is  a  full  three  time  level  scheme. 
n,m  n,m  n,m 


In  order  to  develop  a  three  time  level  volume  consistent  transport  scheme,  we 

;  e.g.  , 


k  k  k 

associate  in  the  advective  terms  d  S  =  d 
-  n  ,m  11,01  ntm 
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38.  The  stability  properties  of  the  above  scheme  were  investigated  for 
the  range  of  conditions  to  be  simulated  in  Mississippi  Sound.  The  scheme  was 
stable  over  this  range  of  flow  conditions.  Details  may  be  found  in  Schmalz 
(1984).  The  coding  of  the  three  time  level  scheme  is  presented  in  Subroutine 
CONCE  in  Appendix  A. 

Dispersion  coefficient  formulation 

39.  To  close  the  numerical  approximations  to  the  two-dimensional,  depth 
averaged  transport  equation,  relations  for  the  effective  dispersion  coeffi¬ 
cients  may  be  developed  in  terms  of  flow  field  properties. 

40.  The  effective  disperison  coefficients  are  assumed  to  have  the  fol¬ 
lowing  form: 


K* 

x 


+  D 

y 


(67) 


where 


K*,K* 
x  y 


u,v 


C 

x 

D 

x 


effective  dispersion  coefficients  in  the  x-  and  y-directions , 
respecitvely 

acceleration  due  to  gravity 

velocity  components  in  the  x-  and  y-directions,  respectively 
water  depth 
Chezy  coefficient 

dispersion  factors  in  the  x-  and  y-directions,  respectively 

dispersion  offsets  due  to  wind  effects  in  the  x-  and  y- 
directions,  respectively  (D  ,D  >  0) 


For  a  unidirectional  flow  in  an  infinitely  wide  channel  in  the  x-direction  , 

Elder  (1959)  found  C  =  5.93  and  C  =  0.23  .  Harleman  et  al.  (1959)  has 

x  y 

converted  Taylor's  result  (1954)  for  pipe  flow  and  determined  =  14.3^2  . 

In  attempting  to  apply  these  results  to  a  two-dimensional  flow  problem  the  fol 

lowing  approach  is  employed.  Initially,  C  ,  C  ,  D  ,  and  D  are  speci- 

x  y  x  y 

fied  by  the  user  as  model  input.  The  cell  face  conditions  for  each  cell  are 

examined  independently  in  each  coordinate  direction.  For  a  no-flux  cell  face 

condition,  C  or  C  and  D  or  D  are  set  to  zero.  For  a  standard  flow 
x  y  x  y 

condition,  the  advective  flag  system  is  examined  to  determine  if  the  flow  is 
restricted  in  the  x-  or  y-direction.  If  the  flow  is  restricted,  C  or 
C  is  reduced  by  a  user-specified  factor. 


PART  III:  MODEL  INPUT  REQUIREMENTS 


41.  The  constituent  transport  schemes  are  included  with  the  hydrody¬ 
namics  as  separate  subroutines  in  WIFM-SAL.  Therefore,  the  model  user  must 
also  be  concerned  with  both  the  hydrodynamic  input  requirements  as  well  as 
those  of  the  transport  computations.  The  complete  input  requirements  for 
WIFM-SAL  are  presented  in  Appendix  B  and  consist  of  29  separate  card  groups. 
Constituent  transport  input  requirements  consist  of  the  following  categories: 

a.  Constituent  Simulation  Control. 

b.  Boundary  Condition  Control. 

c.  Boundary  Condition  Data. 

d.  Wind  Data. 

e.  Constituent  Initial  Condition  Data. 

f.  Dispersion  Coefficient  Data. 

g.  Output  Control. 

Each  category  will  be  discussed  in  detail  below  with  reference  to  the  appro¬ 
priate  card  groups  contained  in  Appendix  B. 

Constituent  Simulation  Control 

42.  This  data  group  is  contained  in  Card  Group  2a.  The  model  user  sets 
ISAL  =1  to  consider  constituent  transport  in  conjunction  with  the  hydrody¬ 
namics.  The  desired  transport  scheme  is  selected  by  specifying  ISALS.  For 
ISALS  =  1  ,  the  FCT  scheme  is  employed,  while  for  ISALS  =  2  ,  the  full  three 
time  level  explicit  scheme  is  used.  Constituent  transport  computations  are 
initiated  ISALC  time  steps  after  the  start  of  the  hydrodynamic  computations. 
The  user  may  set  ISALC  ^  0  in  order  to  allow  for  the  hydrodynamic  computa¬ 
tions  to  be  free  from  initial  condition  effects  before  considering  constituent 
transport.  CMAX  and  XMS  are  self-explanatory. 

Boundary  Condition  Control 

43.  This  data  group  is  contained  in  Card  Groups  3a  and  3b.  In  Card 
Group  3a  the  user  specifies  the  number  of  tidal  elevation  signals  specified  by 
tidal  constituents.  For  a  simulation  over  a  global  grid,  NGLOB  =  0  ,  and  NTI 
is  specified  as  the  number  of  tidal  boundary  (water  surface  elevation  and 


constituent  level)  signals  along  the  seaward  boundary  used  for  interpolation. 
For  a  simulation  over  a  refined  grid,  NTI  =  0  ,  and  NGLOB  is  specified  as  the 
number  of  previously  saved  tidal  signals  (water  surface  elevations  and  con¬ 
stituent  levels)  generated  from  a  global  grid  simulation  to  be  used  for  cell- 
centered  interpolation  along  the  boundary  of  the  refined  grid. 

44.  In  Card  Group  3b,  the  user  specifies  the  grid  indices  for  the  grid 
employed  in  the  current  simulation  where  the  known  tidal  signals  are  available 

Boundary  Condition  Data 

45.  Boundary  condition  format  is  specified  in  Card  Group  3.  The  user 
specifies  ITID  as  the  number  of  entries  in  the  tidal  (elevation  and  constitu¬ 
ents  level)  input  and/or  flow  (discharge  and  constituent  level)  input  data 
tables.  The  number  of  time  steps  between  entries  in  these  tables  is  common 
and  is  specified  as  JTID. 

46.  In  Card  Groups  20c  and  21b,  the  constituent  levels  associated  with 
tidal  and  flow  inputs  are  specified,  respectively. 

Wind  Data 

47.  Detailed  requirements  will  not  be  discussed  here.  Let  it  suffice 
to  say  that  wind  conditions  may  want  to  be  considered  when  simulating  constit¬ 
uent  transport.  The  pertinent  input  variables  requiring  specification  are  as 
follows: 

a.  WA,  THETA  in  Card  Group  4. 

b.  NTABLE  in  Card  Group  5. 

c.  WAT^,  THAT^  in  Card  Group  6  (optional  depending  on  wind 
format) . 

Constituent  Initial  Condition  Data 

48.  In  Card  Group  13a,  the  user  specifies  a  single  format  or  combina¬ 
tion  of  formats  to  be  used  for  specifying  the  constituent  initial  condition. 
IDEPTH  specifies  the  number  of  depth  intervals  used  to  interpolate  based  upon 
depth.  If  IDEPTH  f  0  ,  a  set  of  initial  constituent  levels  TMP^  are  associ¬ 
ated  with  depth  values  D^  j  as  specified  in  Card  Group  13b.  IFIELD 
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specifies  the  number  of  patches  in  which  initial  levels  will  be  specified  on  a 
cell-by-cell  basis.  If  IFIELD  f  0  ,  the  limits  of  patch  and  the  individual 
cell  constituent  levels  are  specified  in  turn  for  each  patch  in  Card  Group  13c 
IZONE  specifies  that  a  number  of  zones  in  which  the  initial  constituent  level 
will  be  a  constant  will  be  assumed.  If  IZONE  f  0  ,  the  number  of  zones,  the 
limits  of  the  zone,  and  the  constant  value  of  initial  constituent  level  for 
the  zone  are  specified  in  Card  Group  13d. 

49.  There  is  considerable  flexibility  in  specifying  initial  constituent 
levels.  Each  format  may  be  used  individually  or  to  override  the  previous 
format.  For  example,  the  user  may  specify  the  initial  conditions  using  depth 
interpolation.  In  selected  areas  of  the  grid  where  detailed  information  is 
available,  the  patch  concept  can  be  used  to  override  the  depth  interpolation. 
In  still  different  areas  of  the  grid,  the  zone  concept  can  be  used  to  specify 
a  uniform  level. 

Dispersion  Coefficient  Data 

50.  Dispersion  factors  and  offsets  due  to  wind  effects  are  specified  in 
turn  for  each  coordinate  direction  in  a  zone  format  as  shown  in  Card  Group  13e 

51.  The  reduction  factor  applied  to  the  dispersion  factors  in  cases  of 
flow  restriction  is  specified  in  Card  Group  17b. 


Output  Control 

52.  Snapshots  of  the  entire  constituent  field  are  printed  after  comple¬ 
tion  of  up  to  32  user-specified  time  steps  during  the  simulation.  Time  step 
completion  data  are  read  in  Card  Group  7  in  the  NPRINT  array. 

53.  Alternatively,  the  user  may  examine  constituent  level  histories  at 
NGAGE  locations  at  NFREQ  time  step  intervals  as  specified  in  Card  Group  5. 

The  NGAGE  locations  are  specified  in  terms  of  the  grid  indicies  in  Card 


PART  IV:  APPLICATION  TO  MISSISSIPPI  SOUND 


54.  Both  the  FCT  and  the  three  time  level  schemes  have  been  applied  to 
the  study  of  salinity  distributions  in  Mississippi  Sound  by  Schmalz  (1984). 

The  schemes  were  exercised  on  a  global  grid  and  also  over  a  local  refined  grid. 
Wind  sensitivity  results  for  both  grid  applications  are  presented  in  turn 
below. 

Global  Grid  Results 


55.  The  horizontal  salinity  distribution  was  simulated  within  Missis¬ 
sippi  Sound  and  adjacent  areas  employing  an  exponentially  stretched  global 
grid  as  shown  in  Figure  3.  This  grid  employs  115  x  59  =  6785  cells.  Maximum 
spatial  resolution  (approximately  3500  ft*)  is  obtained  in  the  passes  into 
Mississippi  Sound.  Depths  within  Mississippi  Sound  are  relatively  shallow 
(10-20  ft),  except  in  the  navigation  channels,  which  are  normally  maintained 
at  30  to  35  ft.  As  a  result,  the  gravity  wave  speed  within  the  Sound  is  less 
than  38  fps,  resulting  in  an  explicit  time  step  limit  of  approximately  100  sec. 
All  simulation  employed  a  360-sec  (6  min)  time  step,  resulting  in  a  maximum 
spatial  Courant  number  of  less  than  4  within  the  Sound. 

56.  Hydrodynamics  and  salinity  conditions  over  the  period  20-24  Sep 
1980  were  simulated.  Water  surface  elevations  along  the  seaward  boundary  were 
obtained  from  a  Gulf  Tide  Model  developed  by  Reid  and  Whitaker  (1981).  Salin¬ 
ity  transect  data  were  available  on  20  and  21  September.  These  values  were 
located  on  the  global  grid  and  two  rectangular  areas  were  set  up  in  which  sa¬ 
linity  values  were  visually  interpolated  from  the  located  transect  values. 
National  Marine  Fisheries  data  were  obtained  for  cruises  No.  106  (Apr  1980) 
and  No.  112  (Nov  1980)  of  the  OREGON  II.  These  data  provided  a  general  under¬ 
standing  of  salinity  patterns  in  the  vicinity  of  the  Mississippi  Delta.  A 
deep  sea  vertically  averaged  value  of  36  ppt  was  employed. 

57.  Initial  conditions  were  assigned  in  a  three  step  process  as  shown 
in  Table  3.  In  step  one,  values  were  assigned  based  on  cell  water  depth.  In 
step  two,  salinity  values  were  specified  within  Mississippi  Sound  based  on 
salinity  transect  data.  In  step  three,  initial  salinity  values  within 


*  A  table  of  factors  for  converting  U.  S.  customary  units  of  measurement  to 
metric  (SI)  is  presented  on  page  3. 
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Lake  Borgne  were  specified  in  a  zone  format.  In  this  process,  each  succeeding 
step  overrides  the  previous  step  values. 

Table  3 

Initial  Salinity  Conditions  on  the  Global  Grid 


Water  Depth 
ft 

Initial  Salinity 
Value,  ppt 

0-10 

22.0 

10-20 

23.0 

20-30 

25.0 

30-50 

30.0 

50-75 

34.0 

75-100 

34.3 

100-120 

34.5 

120-200 

35.0 

200-300 

35.5 

300-500 

36.0 

Salinity  Grid- 

Cell-by-Grid-Cell  Interpolated  Limits 

Global 

Grid  Cell  Range 

Patch 

N 

M 

1 

15-27 

19-39 

2 

28-87 

15-32 

Salinity 

Zone  Specified  Initial  Conditions 

Zone 

Global  Grid 

Cell  Range 

N  M 

Salinity 

ppt 

1 

1-15  33-50 

15 

58.  Salinity  boundary  conditions  which  remained  constant  over  time  are 
shown  in  Table  4.  A  cell-centered  spatial  interpolation  similar  to  that  em¬ 
ployed  for  water  surface  elevations  was  used  to  determine  salinity  values 
along  the  seaward  boundary. 
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Table  4 

Global  Grid  Boundary  Salinity  Conditions 


Tidal  Signal 

Global  Grid  Cell 

Salinity  Value,  ppt 

1 

(115,58) 

36 

2 

(115,56) 

36 

3 

(115,50) 

36 

4 

(115,37) 

36 

5 

(115,22) 

34 

6 

(31,59) 

30 

7 

(42,59) 

36 

8 

(57,59) 

36 

9 

(73,59) 

36 

10 

(87,59) 

36 

11 

(103,59) 

36 

12 

(110,59) 

36 

13 

(112,59) 

36 

14 

(115,59) 

36 

Freshwater 

Inflow 

1 

(97,3) 

0 

2 

(59,19) 

24 

3 

(59,17) 

24 

4 

(13,33) 

15 

5 

(19,20) 

17 

6 

(32,15) 

23 

59.  Wind  data  reported  by  Raytheon  Ocean  Systems  (1981)  over  the  period 
are  presented  in  Table  5.  The  spatially  averaged  wind  speeds  and  directions 
shown  were  lagged  6  hr  in  order  to  investigate  model  sensitivity  to  wind.  A 
constant  drag  coeffecient  equal  to  0.001  was  used  in  the  computations. 

60.  A  total  of  1200  time  steps  were  used  to  simulate  120  hr  of  proto¬ 
type  time.  Wind  information  input  at  6-hr  intervals  was  interpolated  in  time 
at  each  time  step.  Both  salinity  schemes  were  considered.  The  scheme  1  FCT 
results  and  the  scheme  2  three  time  level  results  are  shown  in  Table  6.  The 
following  previously  calibrated  effective  dispersion  coefficients  are  employed 


Reduction  factor  =  0.0388 


Wind  Data  for  20-24  Se 
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MET  3  and  MET  1  are  "land"  stations. 

MET  4  and  MET  5  are  "island"  stations. 

Speed  (MPH).  W  270 

Direction  (Magnetic). 


Table  6 

Global  Grid  Wind  Sensitivity  Simulation 


20/21  Sep  1980 

Transect  Global  Initial 

Station  Grid  Cell  Measured  Condition 


15,39 

16.0 

16. 

16,35 

17.0 

17. 

16,38 

17.3 

17. 

18,33 

17.5 

17. 

18,38 

19.2 

19. 

20,31 

19.2 

19. 

21,35 

23.7 

24. 

23,29 

21.8 

22. 

24,33 

24.9 

25. 

26,29 

22.0 

22. 

27,24 

22.4 

22. 

27,33 

26.8 

27. 

29,26 

23.7 

24. 

29,20 

24.0 

24. 

31,23 

24.8 

25. 

32,26 

25.6 

26. 

33,29 

27.3 

27. 

34,23 

25.2 

25. 

34,31 

28.3 

28. 

34,32; 

28.3 

28. 

40,27 

26.1 

26. 

49,2i: 

23.6 

24. 

49,24 

26.9 

27. 

49,27: 

28.2 

28. 

49,29' 

28.3 

28. 

53,25' 

26.3 

26. 

57,28 

27.3 

27. 

59,21 

27.7 

28. 

60,23' 

28.5 

28. 

62,22 

27.3 

27. 

62,24 

28.1 

28. 

62,28 

29.1 

29. 

62,32 

29.7 

30. 

67,26 

27.9 

28. 

71,28 

28.4 

28. 

75,26 

28.1 

28. 

75,30 

28.7 

29. 

76,25 

26.6 

27. 

81,25 

22.5 

22. 

86,25 

22.9 

23. 

_ 24  Sep  1980 _ 

Computed 

Measured  1  2 


4.2 

19 

7.2 

17 

5.1 

17 

7.6 

17 

9.3 

20 

In  regions  of  the  Sound,  the  scheme  1  and  scheme  2  results  are  nearly  identi¬ 
cal  and  are  in  agreement  with  the  calibration  simulation  and  measured  salinity 
values.  However,  in  the  vicinity  of  the  upper  Mobile  Bay  freshwater  inflow, 
the  results  diverge  as  shown  in  Table  7.  The  scheme  1  results  are  nonnegative 
and  exhibit  no  oscillations.  The  scheme  2  results  exhibit  oscillations  behind 
the  freshwater  front. 

Refined  Grid  Results 

61.  In  order  to  investigate  the  salinity  distribution  in  the  vicinity 
of  the  Pascagoula  Channel,  the  refined  grid  shown  in  Figure  4  was  developed. 
This  grid  employs  49  x  28  =  1372  cells.  Maximum  spatial  resolution  of  300  ft 
is  employed  to  represent  the  navigation  channels.  The  configuration  of  the 
channel  system  is  idealized  in  the  grid  in  order  to  reduce  the  number  of  grid 
cells.  A  60-sec  time  step  was  used,  resulting  in  a  maximum  spatial  Courant 
number  of  less  than  8  within  the  grid  system. 

62.  The  20-24  Sep  1980  period  with  6-hr  lagged  wind  considered  on  the 
global  grid  was  studied  on  the  refined  grid.  The  salinity  values  computed  in 
the  global  grid  scheme  1  FCT  simulation  were  saved  and  interpolated  temporally 
and  spatially  to  provide  the  boundary  conditions  for  the  refined  grid  simula¬ 
tion.  Initial  conditions  over  the  refined  grid  were  determined  from  transect 
data  and  input  cell  by  cell.  Zero  salinity  values  for  the  Pascagoula  River 
System  were  input  for  cells  (8,1)  and  (16,1)  in  order  to  establish  a  fresh¬ 
water  front. 

63.  A  time  step  of  60  sec  was  employed  7200  times  in  order  to  simulate 
120  hr  of  prototype  time.  Wind  information  input  at  6-hr  intervals  from 
Table  5  was  interpolated  in  time  at  each  time  step.  The  scheme  1  FCT  scheme 
was  selected  based  upon  its  superior  performance  on  the  global  grid.  Wind 
lagged  simulation  results  using  the  calibrated  effective  dispersion  coeffi¬ 
cients  in  paragraph  60  are  nearly  identical  to  the  calibration  simulation  re¬ 
sults  and  correspond  to  measured  values  as  shown  in  Table  8.  In  order  to  ob¬ 
tain  an  estimate  of  the  freshwater  influence  and  movement  of  the  front,  the 
salinity  field  at  the  end  of  the  simulation  is  shown  in  Table  9.  Note  the 
scheme  1  results  are  nonnegative  and  exhibit  no  oscillation.  The  flow  pattern 
in  the  vicinity  of  the  freshwater  inflow  at  (16,1)  is  extremely  complex.  The 
averaging  processes  employed  in  coupling  the  two  time  level  scheme  1  FCT  with 


the  three  time  level  hydrodynamics  may  contribute  to  the  unusual  distribution 
over  cells  (15-17,1).  These  effects  are  usually  local  and  the  two  time  level 
scheme  1  FCT  resolves  the  edge  of  the  freshwater  front.  Additional  research 
is  warranted  to  flux-correct  scheme  2  thereby  eliminating  the  above  averaging 
of  hydrodynamic  variables  necessary  in  scheme  1. 

64.  The  input  data  for  this  simulation  are  presented  in  Appendix  C. 
Typical  output  from  the  salinity  computations  embodied  within  WIFM  is  shown  in 
Appendix  D. 


Computer  Requirements 

65.  The  resources  required  for  both  the  hydrodynamics  and  the  salinity 
computations  are  shown  in  Table  10.  Scheme  1  is  more  accurate  but  requires 
nearly  three  times  more  computer  time  than  does  scheme  2.  In  general  a  very 
large  scientific  computation  oriented  machine  should  be  utilized  for  applica¬ 
tions  employing  the  number  of  cells  in  the  Mississippi  Sound  study. 

Table  10 

CRAY  1-S  Requirements 


Simulation 

Grid 

Number  of 
Time  Steps 

Total  Field 
Length  (Octal) 

CPUS 

Scheme  1 

(FCT) 

Global 

1200 

1606064 

619 

Scheme  1 

(FCT) 

Refined 

7200 

776152 

777 

Scheme  2 

(Three  Time  Level) 

Global 

1200 

1601756 

340 

Scheme  2 

(Three  Time  Level) 

Refined 

7200 

776266 

475 

66.  The  job  control  language  (JCL)  for  a  global  grid  and  refined  grid 
simulation  is  presented  in  Appendix  E.  It  should  be  noted  that  the  JCL  shown 
is  for  the  CRAY  I-S  Cray  Operating  System  1.09  as  implemented  at  Kirtland  Air 
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APPENDIX  A:  SUBROUTINE  LISTINGS 


A1 
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TPS>=.5«  iSL 


Cf -AiLi’f’&itiW  ot/?t/<-3  rFT  icff « o*/l«/R3i  face 
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LFjKCUT  ir.E 
AKAhE  T  th  ( 


CFT  10-.M04/14/8J1  F»OC  14 


/ion. 


Cfc/ib/oJ  I|S<  0«:n:«l  U1  10SM04/14/4l»  PAbC  lb 


Card  Group  (Format) 

Variable 

Description 

1 

Required  NDTAP 

Input  tape  unit 

(1615) 

la 

Required  ITL 

Identification  title  card,  up  to 

(8A8) 

64  character,  the  1st  8  are  the 
plot  identification 

2 

Required  NMAX 

Horizontal  grid  dimension  (i.e., 

(1615) 

number  of  cells  in  the  n- 
direction) 

MMAX 

Vertical  grid  dimension,  number  of 
cells  in  the  m-direction 

INITL 

0 — initial  condition 

Restart  conditions: 

system  geometry  and 

1 — restart  conditions  (omit  card 

boundary  imput  tables 

ARE  NOT  read  in 

groups  6,  13-18,  and  23) 

hot  start  conditions: 

INITL  =  0  system 

-m--as  for  0,  but  saves  restart 

geometry  and  boundary 

input  MUST  be  read 

data  every  m  tau 

in;  r|  ,  u  ,  and  v 

are  all  that  has 

been  saved 

I  OVER 

Control  variable 
l--simulation 

0--reads  input  only 

IFLVL 

0--flow  formulation 
l--velocity 

LEVEL 

Number  of  time  levels 

ISURG 

0 — tidal  circulation 

1- -storm  surge-horizontal 

coastline 

2- -storm  surge-vertical  coastline 

IFETR 

0--no  feathering  of  tidal  eleva¬ 
tion,  boundary  elevations,  and 
freshwater  discharges 
l--featheri ng  of  the  above 
quantities 

I  HOT 

0--normal  run 

-l--hot  start  information  previously 
saved  for  r),u,v  ifut  <  oii- 
ditions  on  logical  unit  2  -ill 
be  used 

n--save  hot  start  conditions  at 
itime  =  n  on  logical  unit  21 
Note  r)  surface  elevation 

u  is  x-component  of 
veloci ty 

v  is  y-component  of 
velocity 

itime  is  the  number  of  time 
steps  elapsed 


B3 


Card  Group 

(Format) 

Variable 

Description 

2a 

Required 

ISAL 

l--salinity  simulation 

(3I5.2F10.0) 

0--not  simulating  salinity 

Use  a  blank  card  for  this 

I  SALS 

1--FCT  scheme 

group  when 

you  are  not 

ISALC 

2--3  Time  Level  Explicit  Scheme 

simulating 

salinity 

number  of  time  steps  into  the 
simulation,  when  salinity 
starts 

CMAX 

Maximum  salinity  concentration 
allowable,  in  ppt.  If  CMAX  is 
exceeded,  error  message  is  output 

XMS 

Scale  factor  by  which  salinity 
concentrations  are  multiplied  for 
printout  (dimensionless) 

3 

Required 

ITID 

Number  of  entries  in  tidal  input 

(1615) 

table  or  flow  input  table 

JT1D 

Number  of  t's  between  entries  in 
the  tidal  or  flow  input  tables. 
Note:  if  tidal  constituents  are 

used,  set  ITID  equal  to  the  number 
of  t's  in  the  tidal  scenerio. 

Set  ITJD  =  1  .  Cannot  mix  con¬ 
stituent  or  tabular  entries. 

NTID 

Number  of  distinct  tidal  inputs 
(total  number) 

NFLO 

Number  of  distinct  discharge  input 
(total  number) 

Pertains  only 

^•"NPl 

These  are  print  controls  for  the 

to  velocity 

NP2 

output  grid--grid  is  printed  from 

grid 

^  NP3 

N  =  NP1  to  N  =  NP2  in  steps  of 
NP3 

NP1  and  NP2  are  horizontal  indices 
All  vertical  values  for  each  N 

NP3  is  the  increment 

NPR 

Overrides  NP1,NP2,NP3 

2--print  full  grid  of  r|  only 

-2 — print  full  grid  for  H>u»v 
1 — print  from  NP1  to  NP2,  r|  only 
-l--print  from  NP1  to  NP2,  H>u>v 

MPR  Additional  print  control 

1 —  print  flag  arrays  only 

-1 — print  flag  arrays,  flood, 

barrier,  and  tidal  or  flow  data 

2 —  print  flag  arrays,  depths,  and 
Chezy 

-2--print  all 

MSURF  Counter 

Prints  surface  elevation  and  dis¬ 
charge  in  increments  of  the  values 
of  MSURF 


B4 


Variable 


Card  Group  (Format) 


3 

(1615) 

(Continued) 


3a  Required 

(1615) 


3b  Optional 

(1615) 

NT  =  NT I  +  NGLOB 
Omit  if  NT  =  (J) 


3c  Use  only  if  NCON.GT.(|> 


Omit  if  NCON.EQ.<j> 


3d 

(3711) 

Omit  if  NCON.EQ.<|> 


3e 

(1615) 

Omit  if  NCON.EQ.tp 


Description 


KS1 ,KS2 ,KS3 ,KS4  are  flooding  control 

KSI  m--hold  cell  face  CLOSED  for  mi's 

KS2  m--hold  cell  face  OPEN  for  mt's 

KS3  m--hold  SUBMERGED  barrier  charac¬ 

teristic  for  mi's 

KS4  m--hold  OVERTOPPING  BARRIER  charac¬ 

teristic  for  mi's 

KS5  Leave  blank,  not  used  at  present 

KS6  m--updates  wind  routine  every  mt's 

NCON  Number  of  tide  gages  for  which 

tidal  constituents  must  be 
specified 

NGLOB  Number  of  global  grid  tidal  bound¬ 
ary  signals  used  for  cell-centered 
interpolation  along  the  refined 
grid  boundary 

NT1  Number  of  tidal  boundary  signals 

used  for  cell-centered  interpola¬ 
tion  along  the  seaward  boundary 

Location  expressed  in  m  (grid 
coordinate)  for  each  tidal  signal, 
i  =  1,NT 

Location  expressed  in  m  (grid 
coordinate)  for  each  tidal  signal, 

-i  =  1,NT 


IYEAR 

Start 

time 

of 

simulation 

I MONTH 

Start 

time 

of 

simulation 

I  DAY 

Start 

time 

of 

simulation 

I  HR 

Start 

time 

of 

simulation 

IGX. 

l 

IGY. 

i 


ICONST  Refers  to  array  NCONST  which  con¬ 
tains  37  constituents.  To  choose 
how  many  constituents  you  wish  to 
consider,  code  this  variable 
l--consider 
0--skip 


j  varies  from  1  to  NCONT  where 
NCONT  is  the  element  number  of  the 
specific  constituent  you  want  con¬ 
sidered  from  array  NCONST 


NC. 

J 


Card  Group  (Format) 


Variable 


Description 


Example : 


Required 

See  IXPAN  in  card  group  5.  h 
Code  DX  and  DY .  1  map  inch  / 

=  X  number  of  feet  when  ( 
card  group  6  will  be  \ 

utilized 

Map  scale  1:40000 

CX  =  3333.  ( 


4  Required  TAU  Time  step  length,  i.e.  At  (sec) 

(8F10.0) 

Note:  See  IXPAN  in  card  group  5.  /* DX  Vertical  spatial  stepsize  (minimum 

Code  DX  and  DY.  1  map  inch  /  stepsize  for  Of  space)  from  map 

=  X  number  of  feet  when  (  scale  use  1  in.  =  ft 

card  group  6  will  be  \  „  ,  .  ,  , „  . 

utilized  '  °Y  Horizontal  spatial  stepsize  (ft) 

Example:  Map  scale  1:40000  1  in.  -  -  ft 

CX  =  3333.  G  Acceleration  of  gravity  set  to 

32.2  (ft2/sec) 

ALAT  Average  latitude  of  the  study 

region,  +  for  Northern  hemisphere 
(•)  -  for  Southern  hemisphere 

XI  Constant  rate  of  rainfall  (inches/ 

day) 

WA  Constant  wind  velocity  (no  N/S 

wind) 

-Invariable  wind  as  a  function  of 
time  only.  (Note:  card 
group  11  is  needed  to  complete 
wind  information) 

-2--variable  wind  as  a  function  of 
space  and  time.  (Provided  by 
subroutine  FETCHW.)  (Note: 
omit  card  group  11) 

THETA  Constant  wind  direction  in  degrees 
Use  meteorological  definition,  i.e 
NORTH  is  0° 

EAST  is  90° 

SOUTH  is  180° 

WEST  is  270° 

or  the  number  of  hours  between 
entries  of  wind  table  if  WA  =  -1 

EPSD  is  minimum  amount  of  water 

defining  a  dry  cell  (in  feet) 

APSD  is  minimum  amount  of  water 

over  a  barrier  for  submergence 
(in  feet) 

DC0N1  Value  to  add  to  water  depths  to 
NVGD  translate  them  to  the  model  datum 

1  which  is  usually  NVGD  datum  (in 

MLW  feet).  Depths  are  negative,  thus 

a  -  DC0N1  will  deepen 

*DMPX  Value  of  land  elevation  assigned 
artificially  to  areas  that  will 
never  flood  (in  feet)  control 
value  to  cutoff,  depth  checking 
within  WIFM;  i.e.,  this  is  the  MAX 
land  elevation  digitized 


*  Related  to  XLAND  as  follows:  DMPX  is  used  in  digitizing  the  grid — XLAND  <  DMPX 
defines  highest  potential  flood  level  elevation. 


THETA 


DC0N1 

NVGD 

1 

MLW 

*DMPX 


Card  Group  (Format) 


Variable 


Description 


4 

(8F10.0) 

(Continued) 


Note:  XLAND  <  DMPX  defines 
maximum  potential 
flood  level  elevation 


ROTA  Angle  of  x-axis  as  measured  coun¬ 

terclockwise  from  EAST  =  0°  (in 
degrees) 

TPRO  Start  of  prototype  time  for  be¬ 

ginning  of  run  (i.e.,  time  of  day 
in  hours) 

ADV  0--no  advective  or  viscosity  terms 

1  —  include  advective  terms, 
linearize  at  boundary 
2 — include  advective  terms,  use 
approximation  at  closed  bounds 

VIS  e--viscosity  coefficient  multi¬ 

plier;  it  is  dimensionless  and  if 
equal  to  $  omits  the  viscosity 
coefficient  usually  set  to  1  for 
initial  runs 


XLAND 

NVGD  I  o 

\ 


A  value  of  h  (i.e.,  land  or 
water  bottom  elevation  with  respect 
to  NVGD  datum) ;  greater  than  XLAND 
defines  a  cell  that  will  never 
flood  (in  feet)  XLAND  >  0 


XSCOUR  A  value  of  h  <  XSCOUR  defines  a 

cell  that  will  never  go  dry  (feet) 
0  <  XSCOUR 


SMAX  If  n  >  SMAX  ,  cease  computation 

and  print  q  (H  is  surface  ele¬ 
vation)  (ft) 

SINIT  Set  r|  =  SINIT  as  initial  condi¬ 
tions  (normally  Q) .  Note:  SINIT 
=  999  ,  the  code  will  compute  in¬ 
verted  barometer  effect  (ft) 

DMAXG  Positive  bound  on  maximum  total 
water  depth  that  will  be  ex¬ 
perienced  during  simulation  (in 
feet)  (for  control  of  length  of 
friction  table) 

DC0N2  Value  to  add  to  tidal  input  values 
NVGD  to  translate  them  to  model  datum 
(NVGD)  in  feet 

msl 


DLIMIT  Negative  value  serving  as  an 

artificial  cutoff  value  on  water 
depths  (h)  (negative  since  h  <  0) 
in  feet 


5 

(16IS) 


Required 


MAXTIM  Number  to  t's  to  run  simulation 

INTAP  m — save  q,u,v  on  logical  unit  1 

every  mt 

-l--no  data  is  saved 


Card  Group  (Format) 


5 

(16IS) 

(Continued) 


Note:  Set  these 
variables  to  zero; 
subroutines  to 
accomplish  printer 
plots  have  been 
removed  from  the 
program,  but  can 
be  supplied  upon 
District  request 


Variable 


Description 


IDELAY 


T 

If  these 
plug  con¬ 
trols  are 
set  to  zero, 
omit  card 
group  10 

A 


I  PLOT 

IVPLOT 

ICPLOT 


IXPAN 


NGAGE 


NFREQ 

KREST 

NZP 


Delay  saving  data  on  logical 
unit  1  until  ITIME  =  IDELAY 
(Note:  ITIME  counts  the  number  of 
time  steps) 

^0--printer  plots  of  elevation 
hydrographs  will  be  made 
0 — no  plots 

^0--printer  plots  of  velocity  mag¬ 
nitude  hydrographs  will  be  made 
0--no  plots 

^0--printer  plots  of  peak  surge 
elevation  along  the  coast  will 
be  made 
0--no  plots 

^0--read  in  variable  grid  expan¬ 
sion  coefficients  in  card 
group  6  which  will  be  the 
output  file  from  program  GRID 
saved  on  tape  7 
0--indieates  constant  spatial 
step  input  this  step  size  in 
card  group  4  in  DX  and  DY 
variables 

Number  of  locations  where  you  want 
data  saved,  omit  card  groups  8  and 
9  if  NGAGE  =  0  .  Card  group  8 
is  gage  locations  if  NGAGE  =  0  , 
NFREQ  =  0 

Frequency  to  print  hydrodynamics 
at  gage  points  (every  NFREQ  t's) 

Start  run  at  ITIME  =  KREST. 

Set  to  zero  except  for  restart  run 

Number  of  corrections  to  input 
depth  grid;  omit  card  group  14  if 
NZP  =  0 


NZQ  Number  of  corrections  to  input 

coded  friction  grid;  omit  card 
group  16  if  NZQ  =  0 

MDTAP  Logical  unit  for  depth  and  coded 
friction  input  data  (normally  5) 

NTABLE  Length  of  wind  input  data;  i.e., 
number  of  entries  in  the  table 

IGLOB  n--save  boundary  conditions  (on 

logical  unit  25)  from  Global 
Grid  at  n  points  (cells)  for 
later  use  as  forcing  conditions 
to  an  embedded  grid  (will  need 
card  group  29  to  locate  indices) 
0--no  saving  for  later  use 


B8 


Card  Group  (Format) 


(4G20.11) 


variable 


Optional 


This  group  is  created  by 
program  GRID  on  tape  7  and 
is  omitted  if  IXPAN  =  0 
or  INITL  =  1 


7 

(1615) 


8 

(1615) 

Omit  if  NGAGE  =  0 


Required  NPRINT 


Optional 


NPOT . 
i 


MPOT. 

i 


Optional  IGAGEi 


9 

(1615) 

Omit  if  NGAGE  =  0 


v  =  7  (V  +  V  ,  +  V  +  V  ,  - ) 

4  n,m  n-l,m  n,m+l  n=l,m+l 

u  =  7  (U  +U  _  +  U  .  +  V  _  ,) 

4  n,m  n,m-l  n+l,m  n+l,m-l 


U  =  j  (U  +  U  ) 

Z  n,m  n  ,m- 1 


V  =  i  (V  +  V  ,  ) 

2  n,m  n-l,m 


(n,m) 


IT 


10  0; 
(5I5.2F5.0) 

Omit  all  if:  IPLOT  =  0 

ICPLOT  =  0 
and  IVPLOT  =  0 


Optional 


Description 


Dummy  variable  for  the  first  value 
of  GRID  output 

Expansion  coefficients  for  n- 
direction  (horizontal)  of  the  vari 
able  grid  =  1 ,  NYY  NYY  =  2*NMAX 
(dimensionless) 

Expansion  coefficients  for  verti¬ 
cal  direction  (indirect)  of  the 
variable  grid  i  =  1  , 

NXX  NXX  =  2*MMAX  (dimensionless) 


Time  step  index  to  print  grid  an 
array  of  32  elements  thus  allowing 
up  to  32  printouts  (array  must  be 
filled,  so  two  cards  are  required 
to  satisfy  the  read) 

Horizontal  indices  of  locations 
(i.e.  gage)  where  you  want  data 
saved;  location  is  expressed  in 
terms  of  the  horizontal  dimension 
of  the  grid  (N  values)  i  =  1  , 
NGAGE 

Vertical  indices  of  gage  locations 
(M  values)  i  =  1  ,  NGAGE 

Codes  for  methods  of  computing 
flows  at  gage  points  i  =  1 ,  NGAGE 

1- -u,v 

2- -u,v 

3- -u,v  default  l  v  =  4pt  avg  of 

J  v  at  u) 


4- -u,v 

5-  u 

6- -v 


u  =  4  pt  avg 

u  at  v 


These  groups  are  actually  three 
different  sets  of  variables,  each 
set  associated  with  a  type  of 
printer  plot  to  control  format  of 
plots;  variable  list  and  descrip¬ 
tions  are  not  included 

Subroutines  have  been  removed  from 
the  code,  but  can  be  supplied  upon 
District  request 


Card  Group  (Format) 


Variable 


Description 


10a 

(1615) 


11 

(16F5.2) 

Omit  if  WA  NE-1 
ref  group  4 


12 

(10E8.1) 


Optional  IPLOT — controller  for  elevation 

hydrographs 

IVPLOT — controller  for  velocity 
magnitude  hydrographs 
ICPLOT--controller  for  peak  surge 
elevations  (along  the 
coast)  plot 

■*  Ref:  card  group  5 

Optional  WAT.  Variable  wind  velocity  (mph) 

1  i  =  1  ,  NWAT,  NWAT  =  THETA  (see 

group  4) 

THT.  Corresponding  wind  direction  mea- 

1  sured  from  North  as  THT  (deg) 

i  =  1,  NWAT 

Required  This  card  group  codes  terrain  and  barrier 

characteristics.  Each  variable  in  this  card 
group  has  20  values 

XMAN.  Manning's  coefficient  for  each 

code  i  (i  =  1,20)  used  for  de¬ 
fining  friction  (Note:  value  of 
code  (1)  is  used  for  all  water 
outside  the  computational 
boundaries) . 

This  array  must  be  ordered  in  the 
same  manner  as  the  depth  zones  de¬ 
fined  in  card  group  15.  For 
example--lowest  value  to  highest 
value  of  Manning's  coupled  with 
depth  zones  of  deep  to  shallow  (i 
is  dimensionless) 

ZB.  Barrier  height  for  each  code 

1  i  =  1,20  . 

This  array  is  referenced  by  card 
group  17  variable  INDX  (ft) 

CB.  Chezy  coefficient  to  approximate  a 

A  barrier  of  overtipping  for  each 

code  i  =  1,20 

,  r  ,  1/2  v  r  _  1.49,  .1/6 
(Vg  ,  ft  ,  sec)  Cb  -  —  Ub) 

CO.  Admittance  coefficient  for  over- 

1  -  1/2 

topping  barrier  (^g  ,  ft  /sec) 
usual  range  (3-5)  i  =  1,20 

CAYD.  Recession  coefficient  for  draining 

1  of  flood  cell--keyed  by  friction 

codes  (fraction  of  water  depth  to 
be  allowed  to  drain  within  one 

time  step)  i  =  1,20 


rv- .  ■*. 


Card  Group  (Format) 


Variable 


Description 


12 

(10E8.1) 

(Continued) 


13  Required 

(10F8.0) 

Omit  only  if  INITL.EQ.l 


13a  Optional 

(1615) 


13b  Optional 

(16F5.1) 


13c  Optional 

(1615) 


CD. 

l 


CANPY1 . 
i 


CANPY. 


Admittance  coefficient  for  limit¬ 
ing  movement  of  water  onto  flood 
cells — keyed  by  friction  codes 

(Vg  ,  ft^^/sec)  usual  range  (3-5) 
i  =  1,20 

Canopy  coefficients  for  flooding — 
used  to  increase  Manning's  n 
friction  coefficient  over  heavily 
vegetated  marshes. 

(C  dimensionless)  (C-  is  in  feet) 

nc  =  +  Cle  )  f0r 

d  <  5  ft 

Set  Cj  =  0  ,  and  C^  =  1  for 
nonuse .  i  =  1 , 20 


TMP  Depth  grid  array;  depths  at  center 

of  each  grid  cell.  For  row  M  of 
depths  n  =  1  ,  NMAX,  start  a  new 
card  for  each  M  :  units  of  mea¬ 
sure  (ft)  negative  in  sign 


Include  only  if  ISAL  f  0 

IDEPTH  Number  of  depth  intervals  employed 
to  interpolate  initial  salinity 
condition  based  upon  depth 

IFIELD  Number  of  patches  in  which  initial 
salinity  conditions  will  be  input 
on  a  cell-by-cell  basis 

IZONE  Number  of  zones  in  which  the 

initial  salinity  condition  will  be 
a  constant 

Include  if  IDEPTH  f  0 

TMPn  Salinity  initial  value  array, 

W  N  =  1  ,  IDEPTH  (ppt) 

D(N,1)  Depth  value  array,  N  =  1  , 

IDEPTH  +1  (ft) 


Repeat  IFIELD  times.  Include  only  if 

IFIELD  *  0 

NL  Lower  horizontal  limit  of  patch  i 

(n-coordinate  of  cell) 

NU  Upper  horizontal  limit  of  patch  i 

(n-coordinate  of  cell) 

ML  Lower  vertical  limit  of  patch  i 

(m-coordinate  of  cell) 


Bll 


Variable 


Card  Group  (Format) 


Description 


13c 

(1615) 

(Continued) 


MU  Upper  vertical  limit  of  patch  i 

(m-coordinate  of  cell) 

Repeat  (ML  -  MU)  +  1  times 


(16F5.1) 


CN(N,M),  N  =  NL,NU  Initial  salinity  concentration 

(PPO 


13d 


Required 

Omit  only  if  INITL.EQ.l 


13e 


Card  groups  13d  and  13e  are  ref. 
Subroutine  CONST  to  read  in  five 
different  sets  of  values  for  the 
following  conditions: 

CN — initial  salinity  values  (in 
ppt)  required  only  for 
IZONE  f  0 

CX--dispersion  factor  in  the  X- 
dir  (dimensionless) 

DKXX — dispersion  offset  in  the 
X-dir  (ft  /sec) 

CY — dispersion  factor  in  the  Y- 
dir  (dimensionless) 

DKYY — dispersion  offset  in  the 
Y-dir  (ft  /sec) 


The  variables  for  the  card  group  are: 


(4I5.F10.0) 

NZ 

Number  of  zones  covering  the  grid 
1st  card 

(4I5,F10.0) 

NL 

Lower  horizontal  index  of  zone 

2nd  card 

NU 

Upper  horizontal  index  of  zone 

ML 

Lower  vertical  index  of  zone 

MU 

Upper  vertical  index  of  zone 

R 

Value  of  CN,  CX,  DKXX,  CY,  or  DKYY 
to  be  read  in.  This  is  a  single 
value  for  the  set  of  cells  defined 
N  =  L,u  ,  M  =  L,u  ;  i.e.,  cells 
(N,M)  where  NL  <  N  <  Nu  and 

ML  <  M  <  Mu 

Repeat  the  two  cards  of  this  group  until  all 
initial  condition  variables  above  are 
satisfied. 


14  Optional  N 

(2I5.F5.1) 


Omit  if  NZP.EQ.0  or  INITL.EQ.l  M 

DNM 


Corrections  to  individual  cell 
depths 

Horizontal  index  of  cell 

Vertical  index  of  cell 

Corrected  depth  of  cell  (ft)  nega 
tive  in  sign,  digitized  depth 
without  model  datum  correction, 
usually  reference  is  nautical 
charts,  MLW  or  MLLW  Gulf  Coast 
Datum 


Card  Group  (Format) 


Variable 


15  Required  N  If  N  =  77  ,  the  next  card  begins 

(3512)  the  friction  codes  for  all  cells 

in  the  grid 

Groups  15,  15a,  or  15a  alternate  If  N  ^  77  ,  subroutine  FRICTN  is 

are  required  unless  INITL  =  0  called.  N  is  the  number  of 

depths  used  to  develop  the  fric¬ 
tion  codes  and  the  card  which 
follows  begins  the  definition  of 
depth  ranges 

Within  SUBROUTINE  FRICTN 

Required  DP.  Depths  to  define  ranges  of  depths 

10.1)  which  correspond  to  the  Manning's  n 

In  this  group-the  scheme  ^  the  XMAN  array,  used  to  develop 

operates  to  assign  a  Manning's  the  frlctl0n  codes  (see  group  12) 

n  based  on  depth  i  =  1 ,  ND  (Note:  N  =  ND  and 

Omit  if  INITL. EQ.l  LE.21) 

DA's  are  negative  values  and  must 
be  put  in  the  same  relative  order 
as  the  values  of  Manning's  n  in 
XMAN.  (Note:  deep  to  shallow  if 
lowest  to  highest  is  the  order  of 
Manning's  n)  (ft) 

Required  Use  this  group  if  SUBROUTINE  FRICTN  is  not 

12)  used.  This  alternate  group  15a  is  the  fric¬ 

tion  codes  and  is  related  to  card  group  12 
(XMAN) 

In  this  group  the  scheme  ITIDEM  This  variable  is  read  within  a 

operates  to  allow  you  to  DO  LOOP  where  N  =  1,NMAX  ... 

assign  the  Manning's  n  M  =  1,MMAX 

„ . .  .  ,  tmttt  m  ,  ITIDE  is  a  number  for  each  cell  in 

the  grid  between  1  and  20  cor¬ 
responding  to  the  elements  of  array 
XMAN  which  you  wish  to  assign  to 
each  cell.  The  loops  operate  to 
assign  values  by  columns 

Optional  Used  for  corrections  to  coded 

>)  friction  grid  for  individual  cells 

Omit  if  NZQ.EQ.O  or  INITL. EQ.l  N  Horizontal  cell  index 

M  Vertical  cell  index 

MAN(N,M)  Number  1  to  20  to  correspond  to 

the  elements  of  array  XMAN  desired 

The  code  sets  up  flag  areas  (internal  flags)  for  each  u  and  v  cell  face 
(see  diagram  following)  in  the  grid  based  upon  depth  field  and  advection  code 
(ADV  is  set  equal  to  0,  1,  or  2). 

Card  Group  17  provides  for  establishing  the  codes  for  boundary  flags  and  card 
group  18  provides  for  correction  to  the  internal  flags. 


Required 


Groups  15,  15a,  or  15a  alternate 
are  required  unless  INITL  =  0 


15a 

(8F10. 1) 


Required 


In  this  group--the  scheme 
operates  to  assign  a  Manning's 
n  based  on  depth 

Omit  if  INITL. EQ.l 


15a 

(3512) 


Description 


Required 


In  this  group  the  scheme 
operates  to  allow  you  to 
assign  the  Manning's  n 

Omit  if  INITL. EQ.l 


16  Optional 

(415) 

Omit  if  NZQ.EQ.O  or  INITL. EQ.l 


Card  Group  (Format) 


Variable 


17a  Required 

(312,614) 

Omit  only  if  INITL.EQ.l 


M,X 


Description 


ITYP  Barrier  type  codes 

1 —  exposed  barrier  at  all  times 

2- -overtopping  barrier 

4--submerged  barrier 

8- -tidal  input 

9 —  flow  input 

99 — exit  this  group  of  input, 

leave  remainder  of  card  blank 
UNLESS  you  wish  to  make  cor¬ 
rections  to  the  ICU  and  ICV 
flag  arrays  then  leave  INDX 
and  IDIR  blank  and  set  II  /  0 
and  include  card  group  18.  It 
should  be  set  to  the  number  of 
corrections  to  be  made 

INDX  Value  is  from  1  to  20,  keyed  to 

element  of  array  ZB  in  card  group 
12  to  set  barrier  heights  if  ITYP 
is  set  to  1,  2,  or  4 

For  ITYP  set  to  8,  value  is  from 
1  to  NTID  (NTID  is  the  total  num¬ 
ber  of  tidal  input  signals),  i.e., 
identifies  which  tidal  input 

For  ITYP  set  to  9,  value  is  from 
1  up  to  NFLD  to  identify  which 
flow  input 

IDIR  1 — flow  direction  is  through  u 

cell  Sat ce 

2--flow  direction  is  through  v 
cell  face 


Locator  grid  indices  for  barrier, 
tidal  input,  or  flow  input 

For  a  u  face  feature: 

11  is  the  row  (M) 

12  is  the  beginning  column  (N) 

13  is  the  ending  column  (N) 

For  a  v  face  feature: 

11  is  the  column  (N) 

12  is  the  beginning  row  (M) 

13  is  the  ending  row  (M) 


Card  Group  (Format) 


Variable 


Description 


17a 

(312,614) 

(Continued) 


17b 

(F8.0) 

Include  if  ISAL  f  0 


Optional 


Used  with  tidal  or  flow  input 
only,  otherwise  leave  blank 
0--input  directed  toward  the  right 
or  bottom  of  the  grid 
1 — input  directed  toward  the  left 
or  top  of  the  grid 

When  used  ITYP  =  8  ,  INDX  =  1 
Used  for  tidal  input  only  when  you 
want  to  interpolate  the  values  for 
the  tidal  input  boundary  between 
two  tidal  signals.  15  and  16  cor¬ 
respond  to  the  two  tidal  signal 
numbers;  i.e.,  the  elements  num¬ 
bers  (of  your  2  signals)  in  the 
tidal  signal  arrays  IGX  and  IGY 

Dispersion  coefficient  reduction 
factor  for  flow  restriction 
(dimensionless) 


18 

(415) 

Omit  unless 
ITYP.EQ. 99 .AND. II .NE . 0 


Optional 


See  description  at  ITYP  and 
INDX  codes  above 


19 

(4I5.F5.0) 


Required 


KSHFT 


Correction  codes  to  ICU,  ICV  flag 
arrays.  This  read  statement  is  in 
a  loop  which  will  execute  II  num¬ 
ber  of  times 

Horizontal  index  of  cell 


Vertical  index  of  cell 
'ICU„  M--a  2  digit  code,  n  n  , 

’  where  ^  is  ITYP  and  n 2 
»is.  INDX  of  the  specific  u 
face  of  cell  (N,M) 

_ICV„  M — same  except  v  face  condi- 
’  tion  is  described  by  njn2 

This  card  group  is  a  special 
application: 

NBG--set  equal  to  0 
KSHFT — set  equal  to  1 

0--normal 

-1 — no  tidal  input  or  discharge, 
used  for  storm  surge 

Time  index  unit,  where  the  simula¬ 
tion  begins  in  the  boundary  input 
tables;  i.e.,  time  step  index  for 
beginning  of  input  used  with 
HOTSTART  CONDITIONS 


20 

(4Ab,3F10.0) 

Omit  if  NT10  =  0 
or  NCON  =0  or 
NBG  =  -1 


Optional 
(Tidal  inputs) 


TITLE  Gage  title 

TLON  Longitude  in  degrees  of  gage 

TM  Time  meridian  in  degrees 

HO  Mean  value  referenced  to  model 


Optional 


HM.  . 

J.i 


20a 

(8F10.0) 

Use  only  with  20 


Tidal  amplitude  for  each  of  the 
tidal  constituents 
j  =  1,NC0NT  (NCONT  is  the  element 
of  the  array  NCONST 
which  holds  the  values 
of  the  tidal 
constituents) 

i  =  1,NTID  (number  of  tidal  inputs) 
that  is,  HM  is  the  spe¬ 
cific  tidal  amplitude 
in  feet  of  each  con¬ 
stituent  identified  by 
its  element  number  in 
the  array  NCONST  for 
each  of  the  distinct 
tidal  inputs 


20b  Optional 

(15F5.2) 

These  card  groups  are  tidal 
boundary  INPUT,  thus  omitted 
if  NTID  =  0 


KAPPA.  .  Tidal  phases  of  the  constituents 
^,:L  as  above  (in  degrees) 

SSV.  Tidal  elevation  for  each  time  step 

J  j  =  1,IT  (IT  =  ITID  (card  group 

3)  the  number  of  entries 
in  the  tidal  input  table) 

XKQ  Shift  in  time  step  units 

ALP  Amplitude  multiplication  factor 


Repeat  card  group  20b  (NTID-NCON)  times.  Omit  if  NTID  =  NCON 


NOTE:  If  you  are  simulating  salinity  ISAL.NE.O,  you  will  need  the  following 
additional  groups  repeated  (NTID-NCON)  times. 


20c 

(15F5.2) 

SSV. 

J 

j  =  1 ,IT--specif ic  salinity  value 
for  each  tidal  signal 

21a  Optional 

(15F5.2)  omit  if  NFLO  =  0 

These  groups  are  flow  inputs 

SSV. 

J 

Discharge  (in  cfs)  for  each  time 
step 

j  =  1,IT  where  IT  =  ITID  ,  number 
of  entries  in  flow  input 
table 

XKQ 

Shift  in  time  step  units 

ALP 

Multiplication  factor 

Repeat  this  card 

group  for  each  flow 

input 

If  ISAL.NE.O  you 

will  need  group  21a 

repeated 

for  each  flow  input  as  well 

21b 

(15F5.2) 

Optional 

SSV. 

J 

j  =  1,IT  ,  specific  salinity  value 
for  each  flow  input 
(NFLO) 

22 

Optional 

Specify  if  IFETR  t  0 

(F5.0.I5) 

XLEVEL 

Feather  level  for  tidal  elevation 

NTD 

Number  of  time  steps  for  flow  input 
feathering 

Card  Group  (Format) 


Variable 


Description 


Card  Groups  23  and  24  are  presently  not  being  utilized.  Subroutines  have  been  re- 


moved  which  accomplish 

range  output. 

Supplied  upon  request. 

23 

(1615) 

JNS 

Number  of  ranges  for  computing 
volumetric  discharge  (if  equal  to 

0  put  in  a  blank  card  for  this 
group)  integrated  value 

This  card  group  controls  range 
output 

JT1 

Time  index  marking  beginning  of 
discharge  computation  (time  step) 

Omit  if  INITL  =  1 

JPER 

Period  of  discharge  cycle  in  time 
index  units  (total  length  of 
cycles  in  time  steps) 

Use  a  blank  if  JNS 

=  0 

JDT 

Sampling  time  step  in  time  index 
units 

JMUL 

Number  of  seconds  in  sampling 
period  (l  JDT) 

JDELAY  Delay  print  of  special  gage  data 

until  ITIME  =  JDELAYS  (avoid 
spinup  time  computation) 

24 

(1615) 

Omit  if  JNS  =  0 

Optional 

JDIR . 

i 

Direction  of  flow  in  discharge 
range.  Coded: 

1- -vertical  direction 

2- -horizontal  direction 
where  i  =  i,JNS 

JMN . 

X 

Coordinate  index  of  the  range  line 
i  =  1 , JNS 

JMN1 . ) 

l  Range  line  extends  from  JNS1.  to 

JNS2 .  1 

>  l 

JMN 2 .  I 
^  / 

where  i  =  1,JNS 

25 

(1615) 

Required 

MNPOT 

Number  of  special  gage  points  to 
be  punched  and/or  plotted  for  sur¬ 
face  elevation  data 

MSKP 

Frequency  to  punch  surface  eleva¬ 
tion  data  (i.e.,  every  MSKP  t's) 

MDLY 

Delay  punch  of  surface  elevation 
data  (at  special  gage  points)  until 
ITIME  =  MDLY 

NVELPN  Number  of  special  gage  points  for 

punching  and/or  plotting  velocity 
magnitude 

MVELP 

Frequency  to  punch  velocity  mag¬ 
nitude  data  (every  MVELP  T's)  (x 
is  time  step) 

MVDLY 

Delay  punch  of  velocity  data  until 

ITIME  =  MUDLY 

These  control  variables  are  for  tape  (punch)  or  plot  output  and  are  output  to 
tape  3.  The  plotting  is  automatic.  Plot  file  name  must  be  specified  in  JCL. 


Card  Group  (Format) 


Variable 


Description 


26  Optional 

(1615) 

Omit  if  NNPO7.EQ.0 


27  Optional 

(1615) 

Omit  if  NVELPN 
EQ.0 


28  Optional 

(1615) 

Omit  if  ISURG  =  0 
(see  card  group  2) 


INPOT.  N  indices  of  special  gage  points 
for  surface  elevations  data  (ref. 
card  group  25)  i  =  l.NNPOT  up 
to  a  maximum  of  30  pts 

JMPOT.  M  indices  of  same  (start  a  new 
1  card)  i  =  l.NNPOT 

NVCORD.  N  indices  of  special  gage  points 

for  velocity  magnitude  data  i  =  1 
NVELPN 

MVCORD.  M  indices  of  same  (start  new  card) 
1  i  =  1 .NVELPN 

COAST.  Indices  of  open  coast  cells  where 
i  =  l.NMAX  for  horizontal  coast¬ 
line  and  i  =  l.MMAX  for  vertical 
coastline.  This  information  is 
obtained  from  program  SHORE 


29 


Optional  IU1 .  N  indices  of  cells  where  boundary 

conditions  are  to  be  saved  in 
global  grid,  i  =  l.IGLOB 

M  indices  of  same  i  =  l.IGLOB 
(start  a  new  card) 


15) 

Omit  if  IGLOB  =  0 
(see  card  group  5) 


IU2 . 

l 
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OISCMAftfrt  HOW* 


AMNI1»S 


Global  Grid  Hydrodynamics  and  Salinity 


COESR , T1 7 , 1060 , STCRA . 

ACCOUNT (H485 1 1KV00 , 921C0933-AC0 , ACOCOE , 6343809 ) 

JOB , JN=COESR ,T=120 . 

SWITCH, CARET=+. 

ACQUIRE , DN=$PL , PDN=SFMSHYDDU , RT=0 , MF= I 1 , UQ , ID=C0FHO933 . 
UPDATE , C=DATA,E ,DW=80 , IN . 

DELETE, DN=$PL,NA. 

RELEASE ,DN=$PL. 

ACQUIRE , DN=$PL , PDN=CLMSHYDPU , RT=0 , MF=I 1 , UQ , ID=C0EH0933 . 
UPDATE, C,F, IN. 

DELETE, DN=$PL,NA. 

RELEASE, DN=$ PL. 

AUDIT ,PDN=- , ID=C0EH0933 . 

CFT , I=$CPL . 

REWIND , DN=DATA . 

COPYSBF , I=DATA , 0=$0UT . 

REWIND, DN=DATA. 

ASSIGN, DN=DATA,A=FT05 . 

LDR,MAP,LIB=DISSPLA. 

REWIND ,DN=FT07 : FT08 :FT09 : FT10 : FT1 1 : FT12 . 

COPYD , I=FT07 . 

COPYD , I=FT08 . 

COPYD, I=FT09. 

COPYD, I=FT10. 

COPYD, I=FT11. 

COPYD, I=FT12. 

DISPOSE , DN=FT99 , SDN=COETGS 1PLOT , ID=C0EH0933 , DF-SB , DC=ST , WAIT . 
DISPOSE , DN=FT25 , ID=C0EH0933 , DC=ST , DF=TR , WAIT , + 

TEXT= • USERNO , 868 ,XUMP0M . MASS . STORE  FT25 : TINSH1 . ' . 

DISPOSE ,DN=FT35 , ID=C0EH0933 ,DC=SI ,DF=TR ,  WAIT , + 

TEXT= ' USERNO , 868 , XUMPOM . MASS  STORE  FT35.-TINTS1 . '  . 

EXIT. 

REWIND, DN=FT07 :FT08 :FT09 :FT10 : FT1 1 :FT12. 

COPYD, I=FT07. 

COPYD, I=FT08. 

COPYD, I=FT09. 

COPYD, I=FT10. 

COPYD, I=FT11. 

COPYD, I=FT12. 

DISPOSE ,DN=FT89 , SDN=COETGS 1PLOT , ID=C0EH0933 ,DF=SB , DC=ST , WAIT . 
DISPOSE, DN=FT25 , ID=C0EH0933 ,DC=ST , DF=TR , WAIT , + 

TEXT=’ USERNO, 868, XUMPOM. MASS  STORE  FT25:TINSH1. 

DISPOSE , DN=FT35 , ID*C0EH09  33 , DC=ST , DF=TR , WAIT , + 

TEXT=' USERNO, 868, XUMPOM. MASS  STORE  FT35:TINTS1. 


Refined  Grid  Hydrodynamics  and  Salinity 


COESR , T 1 7 , I 060 , STCRA . 

ACCOUNT (H485 1 1KV00 , 92 1C0933-AC0 , ACOCOE ,6343809 ) 

JOB , JN=COESR , T=820 , CL=C . 

SWITCH, CARET=+. 

ACCESS , DN=FT24 , UQ , NA . 

ACCESS , DN=FT34 , UQ , NA . 

DELETE ,DN=FT24,NA. 

DELETE, DN=FT34,NA. 

RELEASE, DN=FT24. 

RELEASE, DN=FT34. 

ACQUIRE , DN=FT24 , ID=COEH09  33 , RT=0 , DF=TR , UQ , + 

TEXT='USERNO, 868, XUMPOM. MASS. GET  FT24:TINSH. ' . 

ACQUIRE , DN=FT34 , ID=COEH09  3  3 , RT=0 , DF=TR , UQ , + 

TEXT  ='USERNO, 868, XUMPOM. MASS. GET  FT34: TINTS. ’ . 

REWIND , DN=FT24 : FT34 . 

ACCESS ,DN=DATA , ID=COEH0933 ,UQ,NA . 

DELETE, DN=DATA,NA. 

RELEASE, DN=DATA. 

ACQUIRE ,DN=DATA , ID=COEH0933 ,RT=0 ,UQ , + 

TEXT= ’ USERNO , 868 , XUMPOM . MASS . GET  DATA : H485TRSD . ' . 

ACQUIRE , DN=$PL , PDN=CLMSHYDPU , RT=0 , MF=I 1 , UQ , ID=C0EH0933 . 
UPDATE, C,F, IN. 

DELETE, DN=$PL,NA. 

RELEASE, DN=$PL. 

AUDIT , PDN=- , ID=COEH0933 . 

CFT , I=$CPL ,L=0 . 

REWIND, DN=DATA. 

COPYSBF , I=DATA , 0=$0UT . 

REWIND , DN=DATA . 

ASSIGN ,DN=DATA ,A=FT05 . 

LDR ,MAP ,LIB=DISSPLA . 

REWIND , DN=FT0  7 : FT08 : FT09 : FT 10 : FT 1 1 : FT 1 2 . 

C0PYD,I=FT07. 

COPYD , I=FT08 . 

COPYD , I=FT09 . 

COPYD, I=FT10. 

COPYD, I=FT11. 

COPYD, I=FT12. 

DISPOSE ,DN=FT99 ,SDN=COERTSPLOT , ID=COEH0933 ,DF=SB ,DC=ST ,WAIT 
EXIT. 

REWIND , DN=FT07 : FT08 : FT09 : FT 1 0 : FT 1 1 : FT 1 2 . 

COPYD, I=FT07. 

COPYD, I=FT08. 

COPYD, I=FT09. 

COPYD, I=FT10. 

COPYD, I=FT1 1 . 

COPYD, I=FT12. 

DISPOSE , DN=FT99 , SDN=COERTSPLOT , ID=COEH0933 , DF=SB , DC=ST , WAIT 


END 

FILMED 

9-85 

DTIC 


